%%
% 巴特沃斯型
clc;
OmegaP = 2 * pi * 3000;
OmegaS = 2 * pi * 6000;
Rp = 2;
As = 30;
[N, OmegaC] = buttord(OmegaP, OmegaS, Rp, As, 's');
[b , a] = butter(N, OmegaC, 's');
w0 = [OmegaP, OmegaS];
[H, w] = freqs(b, a);
Hx = freqs(b, a, w0);
dbHx = -20 * log10(abs(Hx) / max(abs(H)));
%subplot(2, 2, 1);
plot(w / (2 * pi) / 1000, 20 * log10(abs(H)));xlabel('f(KHz)');ylabel('dB');
set(gca, 'xtickmode', 'manual', 'xtick', [0, 1, 2, 3, 5, 6, 7]);
set(gca, 'ytickmode', 'manual', 'ytick', [-40, -30, -20, -10, 0]);
grid;

%%
%  切贝雪夫 I 型
clc;
OmegaP = 2 * pi * 3000;
OmegaS = 2 * pi * 6000;
Rp = 2;
As = 30;
[N, OmegaC] = cheb1ord(OmegaP, OmegaS, Rp, As, 's');
[b , a] = cheby1(N, Rp, OmegaC, 's');
w0 = [OmegaP, OmegaS];
[H, w] = freqs(b, a);
Hx = freqs(b, a, w0);
dbHx = -20 * log10(abs(Hx) / max(abs(H)));
%subplot(2, 2, 1);
plot(w / (2 * pi) / 1000, 20 * log10(abs(H)));xlabel('f(KHz)');ylabel('dB');
set(gca, 'xtickmode', 'manual', 'xtick', [0, 1, 2, 3, 5, 6, 7]);
set(gca, 'ytickmode', 'manual', 'ytick', [-40, -30, -20, -10, 0]);
grid;

%%
clc;
OmegaP = 2 * pi * 1000;
OmegaS = 2 * pi * 1500;
Rp = 1;As = 15;Fs = 10 * 10^3;
wp = OmegaP / Fs;
ws = OmegaS / Fs;
[N, OmegaC] = buttord(OmegaP, OmegaS, Rp, As, 's');
[b, a] = butter(N, OmegaC, 's');
[bz, az] = impinvar(b, a, Fs);
w0 = [wp, ws];
Hx = freqz(bz, az, w0);
[H, w] = freqz(bz, az);
dbHx = -20 * log10(abs(Hx)/max(abs(H)));
[db, mag, pha, grdm w] = freqz_m(bz, az);
plot(w/pi, db);
xlabel('\omega/\pi');
ylabel('dB');
axis([0, 0.5, -20, 5]);
set(gca, 'xtickmode', 'manual', 'xtick', [0, 0.1, 0.2, 0.3, 0.4, 0.5]);
set(gca, 'ytickmode', 'manual', 'ytick', [-20, -15, -10, -5, -1]);
grid;

%%
clc;
Rp=2;
As=20;
Fs=2000;
wp1=2*pi*200/Fs;
wp2=2*pi*300/Fs;
ws1=2*pi*100/Fs;
ws2=2*pi*400/Fs;
OmegaP1=2*Fs*tan(wp1/2);
OmegaP2=2*Fs*tan(wp2/2);
OmegaS1=2*Fs*tan(ws1/2);
OmegaS2=2*Fs*tan(ws2/2);
OmegaP=[OmegaP1,OmegaP2];
OmegaS=[OmegaS1,OmegaS2];
[N,OmegaC]=buttord(OmegaP,OmegaS,Rp,As,'s');
[b,a]=butter(N,OmegaC,'s');
[bz,az]=bilinear_m(b,a,Fs);
w0=[ws1,wp1,wp2,ws2];
Hx=freqz(bz,az,w0);
[H,w]=freqz(bz,az);
dbHx=-20*log10(abs(Hx)/max(abs(H)));
[db,mag,pha,grd,w]=freqz_m(bz,az);
plot(w/pi,db);
xlabel('\Omega/\pi');
ylabel('dB');
axis([0,0.6,-35,5]);
set(gca, 'xtickmode', 'manual', 'xtick', [0, 0.1, 0.2, 0.3, 0.4, 0.5,0.6]);
set(gca, 'ytickmode', 'manual', 'ytick', [-30, -25, -20, -15, -10, -5, -1]);
grid;

%%
clc;
wp1=0.3;wp2=0.7;
ws1=0.4;ws2=0.6;
wp=[wp1,wp2];ws=[ws1,ws2];
Rp=2;As=30;
[N,wc]=ellipord(wp,ws,Rp,As);
[b,a]=ellip(N,Rp,As,wc);
w0=[wp1*pi,ws1*pi,ws2*pi,wp2*pi];
Hx=freqz(b,a,w0);
[H,w]=freqz(b,a);
dbHx=-20*log10(abs(Hx)/max(abs(H)));
[db,mag,pha,grd,w]=freqz_m(b,a);
plot(w/pi,db);title('椭圆函数型数字带阻幅度响应（dB）');
xlabel('\omega/\pi');ylabel('dB');
axis([0,1,-60,4]);grid;

%%
% 设计一个数字巴特沃思高通滤波器
% 衰减As>=30dB,当f<=3kHz
% 衰减Rp<=3dB,当f>=5kHz
% 抽样频率fs=20kHz
% 试用双线性变换法进行设计，最后写出H(z)的表达式并画出其幅频响应特性（dB）
clc;
clear;

% 规格参数
As = 30;
Rp = 3;
fs = 20000;
fp = 5000;
fs1 = 3000;

wp = 2 * pi * fp / fs;
ws = 2 * pi * fs1 / fs;
OmegaP = 2 * fs * tan(wp / 2);
OmegaS = 2 * fs * tan(ws / 2);
[N, OmegaC] = buttord(OmegaP, OmegaS, Rp, As, 's');
[b, a] = butter(N, OmegaC, 'high', 's');
[bz, az] = bilinear(b, a, fs);
[H, w] = freqz(bz, az, 1024, fs);
db = 20 * log10(abs(H));

% 画出幅频响应特性
figure;
plot(w, db);
title('数字巴特沃斯高通滤波器的幅频响应');
xlabel('频率 (Hz)');
ylabel('幅度 (dB)');
grid on;
% 输出H(z)的表达式
Hz = tf(bz, az, 1/fs);
disp('H(z)的表达式:');
Hz

%%
% 设计一个数字切贝雪夫Ⅰ型带阻滤波器
% 衰减As>=30dB,当1kHz<=f<=2kHz
% 波纹Rp<=3dB,当f>=500Hz，f>=3kHz
% 抽样频率fs=10kHz
% 试用双线性变换法进行设计，最后写出H(z)的表达式并画出其幅频响应特性（dB）
clc;
clear;

% 规格参数
As = 30;
Rp = 3; 
fs = 10000;
fp1 = 500; 
fp2 = 3000;
fs1 = 1000; 
fs2 = 2000; 

wp1 = 2 * pi * fp1 / fs;
wp2 = 2 * pi * fp2 / fs;
ws1 = 2 * pi * fs1 / fs;
ws2 = 2 * pi * fs2 / fs;
OmegaP1 = 2 * fs * tan(wp1 / 2);
OmegaP2 = 2 * fs * tan(wp2 / 2);
OmegaS1 = 2 * fs * tan(ws1 / 2);
OmegaS2 = 2 * fs * tan(ws2 / 2);
[N, OmegaC] = cheb1ord([OmegaP1, OmegaP2], [OmegaS1, OmegaS2], Rp, As, 's');
[b, a] = cheby1(N, Rp, OmegaC, 'stop', 's');
[bz, az] = bilinear(b, a, fs);
[H, w] = freqz(bz, az, 1024, fs);
db = 20 * log10(abs(H));

% 画出幅频响应特性
figure;
plot(w, db);
title('数字切比雪夫Ⅰ型带阻滤波器的幅频响应');
xlabel('频率 (Hz)');
ylabel('幅度 (dB)');
grid on;
% 输出H(z)的表达式
Hz = tf(bz, az, 1/fs);
disp('H(z)的表达式:');
Hz

%%
% 设计一个数字带通滤波器
% 抽样频率fs=25kHz，通带截止频率fp1=5kHz，fp2=7kHz，通带衰减Rp=0.5dB，阻带截止频率fu1=3.5kHz，fu2=8.5kHz，阻带衰减As=45dB
clc;
clear;

% 规格参数
fs = 25000; 
fp1 = 5000; 
fp2 = 7000; 
Rp = 0.5;   
fu1 = 3500; 
fu2 = 8500; 
As = 45; 
   
wp1 = 2 * fp1 / fs;
wp2 = 2 * fp2 / fs;
ws1 = 2 * fu1 / fs;
ws2 = 2 * fu2 / fs;
[N, Wn] = cheb1ord([wp1 wp2], [ws1 ws2], Rp, As);
[b, a] = cheby1(N, Rp, Wn, 'bandpass');
[H, w] = freqz(b, a, 1024, fs);
db = 20 * log10(abs(H));

% 画出幅频响应特性
figure;
plot(w, db);
title('数字切比雪夫Ⅰ型带通滤波器的幅频响应');
xlabel('频率 (Hz)');
ylabel('幅度 (dB)');
grid on;
% 输出H(z)的表达式
Hz = tf(b, a, 1/fs);
disp('H(z)的表达式:');
Hz
